Warning in mcse.multi(x, method = method, size = size): Estimated matrix not
positive definite. The chain might be highly correlated or very high
dimensional. Consider increasing the sample size. Using the default batch means
estimator as a substitute.
Warning in mcse.multi(x, method = method, size = size): Estimated matrix not
positive definite. The chain might be highly correlated or very high
dimensional. Consider increasing the sample size. Using the default batch means
estimator as a substitute.
#all three chains combined all converged on the same parametersknitr::include_graphics("/Users/devonderaad/Desktop/ESS.r2.all.png")
Code
#check Gelman-Rubin diagnostic. If PSRF scores are near 1, then we are confident the chain has converged#rep1.r1x<-read.table("~/Desktop/nmel.ceyx.rad/snapp/rep2.log", header=T)xx<-as.matrix(x[451:2001,2:7])stable.GR(xx)
Warning in mcse.multi(x, method = method, size = size): Estimated matrix not
positive definite. The chain might be highly correlated or very high
dimensional. Consider increasing the sample size. Using the default batch means
estimator as a substitute.
#all three chains combined all converged on the same parametersknitr::include_graphics("/Users/devonderaad/Desktop/ESS.r3.all.png")
Code
#check Gelman-Rubin diagnostic. If PSRF scores are near 1, then we are confident the chain has converged#rep3.r1x<-read.table("~/Desktop/nmel.ceyx.rad/snapp/rep3.log", header=T)xx<-as.matrix(x[452:2001,2:7])stable.GR(xx)
Warning in mcse.multi(x, method = method, size = size): Estimated matrix not
positive definite. The chain might be highly correlated or very high
dimensional. Consider increasing the sample size. Using the default batch means
estimator as a substitute.
Warning in mcse.multi(x, method = method, size = size): Estimated matrix not
positive definite. The chain might be highly correlated or very high
dimensional. Consider increasing the sample size. Using the default batch means
estimator as a substitute.
4 Check whether early burst is found in each replicate: Rep1(r1)
Code
#phytools approachlibrary(phytools)
Loading required package: ape
Loading required package: maps
Code
#read in the consensus treetree <- ape::read.nexus(file ="~/Desktop/nmel.ceyx.rad/snapp/rep1.con.tree")#remove root so that only the diversification rate within the northern Melanesian clade is assessedtree<-drop.tip(tree, "margarethae")#set h = the height of the empirical treeh<-max(nodeHeights(tree))h
[1] 3.792403
Code
#set x as 1000 evenly spaced numbers between zero and the height of the input treex<-seq(0,h,by=h/1000)#set birth rate equal to a log-linear rate, following a pure birth predictionb<-(log(Ntip(tree))-log(2))/h#simulate 1000 trees with the same height and tip number as our empirical tree, plus a log-linear diversification ratetrees<-pbtree(b=b, n=Ntip(tree), t=h,nsim=1000 ,method="direct", quiet=TRUE)#calculate 95% CI for the 1000 trees simmed above under a pure birth/death modelobject<-ltt95(trees,log=TRUE,plot=FALSE, mode ="mean")#set up the plotpar(bty="l")#plot background 95% CI based on simulated dataplot(object, labels=FALSE, shaded=TRUE, xaxis ="negative", bg="grey90")axis(side =1, at=seq(-4,0, by=.5))#cover the mean line with a grey linelines(x=object[,1]-max(object[,1]),y=object[,3], type="l", lwd=2, col="grey60")#plot empirical LTTobj2<-ltt(tree, plot=FALSE)#flip x axis values to negativeobj2$times<-obj2$times-max(obj2$times)#plot on top of existing plotplot(obj2, add=TRUE, log.lineages=FALSE, col="black", lwd=4, lty=1)
5 Check whether early burst is found in each replicate: Rep1(r2)
Code
#read in the consensus treetree <- ape::read.nexus(file ="~/Desktop/nmel.ceyx.rad/snapp/rep1.r2/rep1.con.tree")#remove root so that only the diversification rate within the northern Melanesian clade is assessedtree<-drop.tip(tree, "margarethae")#set h = the height of the empirical treeh<-max(nodeHeights(tree))h
[1] 3.788046
Code
#set x as 1000 evenly spaced numbers between zero and the height of the input treex<-seq(0,h,by=h/1000)#set birth rate equal to a log-linear rate, following a pure birth predictionb<-(log(Ntip(tree))-log(2))/h#simulate 1000 trees with the same height and tip number as our empirical tree, plus a log-linear diversification ratetrees<-pbtree(b=b, n=Ntip(tree), t=h,nsim=1000 ,method="direct", quiet=TRUE)#calculate 95% CI for the 1000 trees simmed above under a pure birth/death modelobject<-ltt95(trees,log=TRUE,plot=FALSE, mode ="mean")#set up the plotpar(bty="l")#plot background 95% CI based on simulated dataplot(object, labels=FALSE, shaded=TRUE, xaxis ="negative", bg="grey90")axis(side =1, at=seq(-4,0, by=.5))#cover the mean line with a grey linelines(x=object[,1]-max(object[,1]),y=object[,3], type="l", lwd=2, col="grey60")#plot empirical LTTobj2<-ltt(tree, plot=FALSE)#flip x axis values to negativeobj2$times<-obj2$times-max(obj2$times)#plot on top of existing plotplot(obj2, add=TRUE, log.lineages=FALSE, col="black", lwd=4, lty=1)
6 Check whether early burst is found in each replicate: Rep1(r3)
Code
#read in the consensus treetree <- ape::read.nexus(file ="~/Desktop/nmel.ceyx.rad/snapp/rep1.r3/rep1.con.tree")#remove root so that only the diversification rate within the northern Melanesian clade is assessedtree<-drop.tip(tree, "margarethae")#set h = the height of the empirical treeh<-max(nodeHeights(tree))h
[1] 3.839191
Code
#set x as 1000 evenly spaced numbers between zero and the height of the input treex<-seq(0,h,by=h/1000)#set birth rate equal to a log-linear rate, following a pure birth predictionb<-(log(Ntip(tree))-log(2))/h#simulate 1000 trees with the same height and tip number as our empirical tree, plus a log-linear diversification ratetrees<-pbtree(b=b, n=Ntip(tree), t=h,nsim=1000 ,method="direct", quiet=TRUE)#calculate 95% CI for the 1000 trees simmed above under a pure birth/death modelobject<-ltt95(trees,log=TRUE,plot=FALSE, mode ="mean")#set up the plotpar(bty="l")#plot background 95% CI based on simulated dataplot(object, labels=FALSE, shaded=TRUE, xaxis ="negative", bg="grey90")axis(side =1, at=seq(-4,0, by=.5))#cover the mean line with a grey linelines(x=object[,1]-max(object[,1]),y=object[,3], type="l", lwd=2, col="grey60")#plot empirical LTTobj2<-ltt(tree, plot=FALSE)#flip x axis values to negativeobj2$times<-obj2$times-max(obj2$times)#plot on top of existing plotplot(obj2, add=TRUE, log.lineages=FALSE, col="black", lwd=4, lty=1)
7 Check whether early burst is found in each replicate: Rep2(r1)
Code
#read in the consensus treetree <- ape::read.nexus(file ="~/Desktop/nmel.ceyx.rad/snapp/rep2.con.tree")#remove root so that only the diversification rate within the northern Melanesian clade is assessedtree<-drop.tip(tree, "margarethae")#set h = the height of the empirical treeh<-max(nodeHeights(tree))h
[1] 3.806352
Code
#set x as 1000 evenly spaced numbers between zero and the height of the input treex<-seq(0,h,by=h/1000)#set birth rate equal to a log-linear rate, following a pure birth predictionb<-(log(Ntip(tree))-log(2))/h#simulate 1000 trees with the same height and tip number as our empirical tree, plus a log-linear diversification ratetrees<-pbtree(b=b, n=Ntip(tree), t=h,nsim=1000 ,method="direct", quiet=TRUE)#calculate 95% CI for the 1000 trees simmed above under a pure birth/death modelobject<-ltt95(trees,log=TRUE,plot=FALSE, mode ="mean")#set up the plotpar(bty="l")#plot background 95% CI based on simulated dataplot(object, labels=FALSE, shaded=TRUE, xaxis ="negative", bg="grey90")axis(side =1, at=seq(-4,0, by=.5))#cover the mean line with a grey linelines(x=object[,1]-max(object[,1]),y=object[,3], type="l", lwd=2, col="grey60")#plot empirical LTTobj2<-ltt(tree, plot=FALSE)#flip x axis values to negativeobj2$times<-obj2$times-max(obj2$times)#plot on top of existing plotplot(obj2, add=TRUE, log.lineages=FALSE, col="black", lwd=4, lty=1)
8 Check whether early burst is found in each replicate: Rep2(r2)
Code
#read in the consensus treetree <- ape::read.nexus(file ="~/Desktop/nmel.ceyx.rad/snapp/rep2.r2/rep2.con.tree")#remove root so that only the diversification rate within the northern Melanesian clade is assessedtree<-drop.tip(tree, "margarethae")#set h = the height of the empirical treeh<-max(nodeHeights(tree))h
[1] 3.839811
Code
#set x as 1000 evenly spaced numbers between zero and the height of the input treex<-seq(0,h,by=h/1000)#set birth rate equal to a log-linear rate, following a pure birth predictionb<-(log(Ntip(tree))-log(2))/h#simulate 1000 trees with the same height and tip number as our empirical tree, plus a log-linear diversification ratetrees<-pbtree(b=b, n=Ntip(tree), t=h,nsim=1000 ,method="direct", quiet=TRUE)#calculate 95% CI for the 1000 trees simmed above under a pure birth/death modelobject<-ltt95(trees,log=TRUE,plot=FALSE, mode ="mean")#set up the plotpar(bty="l")#plot background 95% CI based on simulated dataplot(object, labels=FALSE, shaded=TRUE, xaxis ="negative", bg="grey90")axis(side =1, at=seq(-4,0, by=.5))#cover the mean line with a grey linelines(x=object[,1]-max(object[,1]),y=object[,3], type="l", lwd=2, col="grey60")#plot empirical LTTobj2<-ltt(tree, plot=FALSE)#flip x axis values to negativeobj2$times<-obj2$times-max(obj2$times)#plot on top of existing plotplot(obj2, add=TRUE, log.lineages=FALSE, col="black", lwd=4, lty=1)
9 Check whether early burst is found in each replicate: Rep2(r3)
Code
#read in the consensus treetree <- ape::read.nexus(file ="~/Desktop/nmel.ceyx.rad/snapp/rep2.r3/rep2.con.tree")#remove root so that only the diversification rate within the northern Melanesian clade is assessedtree<-drop.tip(tree, "margarethae")#set h = the height of the empirical treeh<-max(nodeHeights(tree))h
[1] 3.798614
Code
#set x as 1000 evenly spaced numbers between zero and the height of the input treex<-seq(0,h,by=h/1000)#set birth rate equal to a log-linear rate, following a pure birth predictionb<-(log(Ntip(tree))-log(2))/h#simulate 1000 trees with the same height and tip number as our empirical tree, plus a log-linear diversification ratetrees<-pbtree(b=b, n=Ntip(tree), t=h,nsim=1000 ,method="direct", quiet=TRUE)#calculate 95% CI for the 1000 trees simmed above under a pure birth/death modelobject<-ltt95(trees,log=TRUE,plot=FALSE, mode ="mean")#set up the plotpar(bty="l")#plot background 95% CI based on simulated dataplot(object, labels=FALSE, shaded=TRUE, xaxis ="negative", bg="grey90")axis(side =1, at=seq(-4,0, by=.5))#cover the mean line with a grey linelines(x=object[,1]-max(object[,1]),y=object[,3], type="l", lwd=2, col="grey60")#plot empirical LTTobj2<-ltt(tree, plot=FALSE)#flip x axis values to negativeobj2$times<-obj2$times-max(obj2$times)#plot on top of existing plotplot(obj2, add=TRUE, log.lineages=FALSE, col="black", lwd=4, lty=1)
10 Check whether early burst is found in each replicate: Rep3(r1)
Code
#read in the consensus treetree <- ape::read.nexus(file ="~/Desktop/nmel.ceyx.rad/snapp/rep3.consensus.tree")#remove root so that only the diversification rate within the northern Melanesian clade is assessedtree<-drop.tip(tree, "margarethae")#set h = the height of the empirical treeh<-max(nodeHeights(tree))h
[1] 3.798202
Code
#set x as 1000 evenly spaced numbers between zero and the height of the input treex<-seq(0,h,by=h/1000)#set birth rate equal to a log-linear rate, following a pure birth predictionb<-(log(Ntip(tree))-log(2))/h#simulate 1000 trees with the same height and tip number as our empirical tree, plus a log-linear diversification ratetrees<-pbtree(b=b, n=Ntip(tree), t=h,nsim=1000 ,method="direct", quiet=TRUE)#calculate 95% CI for the 1000 trees simmed above under a pure birth/death modelobject<-ltt95(trees,log=TRUE,plot=FALSE, mode ="mean")#set up the plotpar(bty="l")#plot background 95% CI based on simulated dataplot(object, labels=FALSE, shaded=TRUE, xaxis ="negative", bg="grey90")axis(side =1, at=seq(-4,0, by=.5))#cover the mean line with a grey linelines(x=object[,1]-max(object[,1]),y=object[,3], type="l", lwd=2, col="grey60")#plot empirical LTTobj2<-ltt(tree, plot=FALSE)#flip x axis values to negativeobj2$times<-obj2$times-max(obj2$times)#plot on top of existing plotplot(obj2, add=TRUE, log.lineages=FALSE, col="black", lwd=4, lty=1)
11 Check whether early burst is found in each replicate: Rep3(r2)
Code
#read in the consensus treetree <- ape::read.nexus(file ="~/Desktop/nmel.ceyx.rad/snapp/rep3.r2/rep3.con.tree")#remove root so that only the diversification rate within the northern Melanesian clade is assessedtree<-drop.tip(tree, "margarethae")#set h = the height of the empirical treeh<-max(nodeHeights(tree))h
[1] 3.82212
Code
#set x as 1000 evenly spaced numbers between zero and the height of the input treex<-seq(0,h,by=h/1000)#set birth rate equal to a log-linear rate, following a pure birth predictionb<-(log(Ntip(tree))-log(2))/h#simulate 1000 trees with the same height and tip number as our empirical tree, plus a log-linear diversification ratetrees<-pbtree(b=b, n=Ntip(tree), t=h,nsim=1000 ,method="direct", quiet=TRUE)#calculate 95% CI for the 1000 trees simmed above under a pure birth/death modelobject<-ltt95(trees,log=TRUE,plot=FALSE, mode ="mean")#set up the plotpar(bty="l")#plot background 95% CI based on simulated dataplot(object, labels=FALSE, shaded=TRUE, xaxis ="negative", bg="grey90")axis(side =1, at=seq(-4,0, by=.5))#cover the mean line with a grey linelines(x=object[,1]-max(object[,1]),y=object[,3], type="l", lwd=2, col="grey60")#plot empirical LTTobj2<-ltt(tree, plot=FALSE)#flip x axis values to negativeobj2$times<-obj2$times-max(obj2$times)#plot on top of existing plotplot(obj2, add=TRUE, log.lineages=FALSE, col="black", lwd=4, lty=1)
12 Check whether early burst is found in each replicate: Rep3(r3)
Code
#read in the consensus treetree <- ape::read.nexus(file ="~/Desktop/nmel.ceyx.rad/snapp/rep3.r3/rep3.con.tree")#remove root so that only the diversification rate within the northern Melanesian clade is assessedtree<-drop.tip(tree, "margarethae")#set h = the height of the empirical treeh<-max(nodeHeights(tree))h
[1] 3.816228
Code
#set x as 1000 evenly spaced numbers between zero and the height of the input treex<-seq(0,h,by=h/1000)#set birth rate equal to a log-linear rate, following a pure birth predictionb<-(log(Ntip(tree))-log(2))/h#simulate 1000 trees with the same height and tip number as our empirical tree, plus a log-linear diversification ratetrees<-pbtree(b=b, n=Ntip(tree), t=h,nsim=1000 ,method="direct", quiet=TRUE)#calculate 95% CI for the 1000 trees simmed above under a pure birth/death modelobject<-ltt95(trees,log=TRUE,plot=FALSE, mode ="mean")#set up the plotpar(bty="l")#plot background 95% CI based on simulated dataplot(object, labels=FALSE, shaded=TRUE, xaxis ="negative", bg="grey90")axis(side =1, at=seq(-4,0, by=.5))#cover the mean line with a grey linelines(x=object[,1]-max(object[,1]),y=object[,3], type="l", lwd=2, col="grey60")#plot empirical LTTobj2<-ltt(tree, plot=FALSE)#flip x axis values to negativeobj2$times<-obj2$times-max(obj2$times)#plot on top of existing plotplot(obj2, add=TRUE, log.lineages=FALSE, col="black", lwd=4, lty=1)
Source Code
---title: "Check that we consistently recover the early burst pattern across replicates and downsampling schemes"format: html: code-fold: show code-tools: truetoc: truetoc-title: Document Contentsnumber-sections: trueembed-resources: true---First, check that the three replicates for each downsampling scheme converged, with appropriate effective sample sizes (>200).## Rep1 check convergence```{r}#show ESS for individual replicates#rep1knitr::include_graphics("/Users/devonderaad/Desktop/ESS.r1.r1.png")#rep2knitr::include_graphics("/Users/devonderaad/Desktop/ESS.r1.r2.png")#rep3knitr::include_graphics("/Users/devonderaad/Desktop/ESS.r1.r3.png")#all three chains combined all converged on the same parametersknitr::include_graphics("/Users/devonderaad/Desktop/ESS.r1.all.png")#check Gelman-Rubin diagnostic. If PSRF scores are near 1, then we are confident the chain has convergedlibrary(stableGR)#rep1.r1x<-read.table("~/Desktop/nmel.ceyx.rad/snapp/rep1.log", header=T)xx<-as.matrix(x[454:2001,2:7])stable.GR(xx)#rep1.r2y<-read.table("~/Desktop/nmel.ceyx.rad/snapp/rep1.r2/rep1.log", header=T)yy<-as.matrix(y[451:2001,2:7])stable.GR(yy)#rep1.r3z<-read.table("~/Desktop/nmel.ceyx.rad/snapp/rep1.r3/rep1.log", header=T)zz<-as.matrix(z[452:2001,2:7])stable.GR(zz)```## Rep2 check convergence```{r}#show ESS for individual replicates#rep1knitr::include_graphics("/Users/devonderaad/Desktop/ESS.r2.r1.png")#rep2knitr::include_graphics("/Users/devonderaad/Desktop/ESS.r2.r2.png")#rep3knitr::include_graphics("/Users/devonderaad/Desktop/ESS.r2.r3.png")#all three chains combined all converged on the same parametersknitr::include_graphics("/Users/devonderaad/Desktop/ESS.r2.all.png")#check Gelman-Rubin diagnostic. If PSRF scores are near 1, then we are confident the chain has converged#rep1.r1x<-read.table("~/Desktop/nmel.ceyx.rad/snapp/rep2.log", header=T)xx<-as.matrix(x[451:2001,2:7])stable.GR(xx)#rep1.r2y<-read.table("~/Desktop/nmel.ceyx.rad/snapp/rep2.r2/rep2.log", header=T)yy<-as.matrix(y[452:2001,2:7])stable.GR(yy)#rep1.r3z<-read.table("~/Desktop/nmel.ceyx.rad/snapp/rep2.r3/rep2.log", header=T)zz<-as.matrix(z[453:2001,2:7])stable.GR(zz)```## Rep3 check convergence```{r}#show ESS for individual replicates#rep1knitr::include_graphics("/Users/devonderaad/Desktop/ESS.r3.r1.png")#rep2knitr::include_graphics("/Users/devonderaad/Desktop/ESS.r3.r2.png")#rep3knitr::include_graphics("/Users/devonderaad/Desktop/ESS.r3.r3.png")#all three chains combined all converged on the same parametersknitr::include_graphics("/Users/devonderaad/Desktop/ESS.r3.all.png")#check Gelman-Rubin diagnostic. If PSRF scores are near 1, then we are confident the chain has converged#rep3.r1x<-read.table("~/Desktop/nmel.ceyx.rad/snapp/rep3.log", header=T)xx<-as.matrix(x[452:2001,2:7])stable.GR(xx)#rep3.r2y<-read.table("~/Desktop/nmel.ceyx.rad/snapp/rep3.r2/rep3.log", header=T)yy<-as.matrix(y[452:2001,2:7])stable.GR(yy)#rep3.r3z<-read.table("~/Desktop/nmel.ceyx.rad/snapp/rep3.r3/rep3.log", header=T)zz<-as.matrix(z[452:2001,2:7])stable.GR(zz)```## Check whether early burst is found in each replicate: Rep1(r1)```{r}#phytools approachlibrary(phytools)#read in the consensus treetree <- ape::read.nexus(file ="~/Desktop/nmel.ceyx.rad/snapp/rep1.con.tree")#remove root so that only the diversification rate within the northern Melanesian clade is assessedtree<-drop.tip(tree, "margarethae")#set h = the height of the empirical treeh<-max(nodeHeights(tree))h#set x as 1000 evenly spaced numbers between zero and the height of the input treex<-seq(0,h,by=h/1000)#set birth rate equal to a log-linear rate, following a pure birth predictionb<-(log(Ntip(tree))-log(2))/h#simulate 1000 trees with the same height and tip number as our empirical tree, plus a log-linear diversification ratetrees<-pbtree(b=b, n=Ntip(tree), t=h,nsim=1000 ,method="direct", quiet=TRUE)#calculate 95% CI for the 1000 trees simmed above under a pure birth/death modelobject<-ltt95(trees,log=TRUE,plot=FALSE, mode ="mean")#set up the plotpar(bty="l")#plot background 95% CI based on simulated dataplot(object, labels=FALSE, shaded=TRUE, xaxis ="negative", bg="grey90")axis(side =1, at=seq(-4,0, by=.5))#cover the mean line with a grey linelines(x=object[,1]-max(object[,1]),y=object[,3], type="l", lwd=2, col="grey60")#plot empirical LTTobj2<-ltt(tree, plot=FALSE)#flip x axis values to negativeobj2$times<-obj2$times-max(obj2$times)#plot on top of existing plotplot(obj2, add=TRUE, log.lineages=FALSE, col="black", lwd=4, lty=1)```## Check whether early burst is found in each replicate: Rep1(r2)```{r}#read in the consensus treetree <- ape::read.nexus(file ="~/Desktop/nmel.ceyx.rad/snapp/rep1.r2/rep1.con.tree")#remove root so that only the diversification rate within the northern Melanesian clade is assessedtree<-drop.tip(tree, "margarethae")#set h = the height of the empirical treeh<-max(nodeHeights(tree))h#set x as 1000 evenly spaced numbers between zero and the height of the input treex<-seq(0,h,by=h/1000)#set birth rate equal to a log-linear rate, following a pure birth predictionb<-(log(Ntip(tree))-log(2))/h#simulate 1000 trees with the same height and tip number as our empirical tree, plus a log-linear diversification ratetrees<-pbtree(b=b, n=Ntip(tree), t=h,nsim=1000 ,method="direct", quiet=TRUE)#calculate 95% CI for the 1000 trees simmed above under a pure birth/death modelobject<-ltt95(trees,log=TRUE,plot=FALSE, mode ="mean")#set up the plotpar(bty="l")#plot background 95% CI based on simulated dataplot(object, labels=FALSE, shaded=TRUE, xaxis ="negative", bg="grey90")axis(side =1, at=seq(-4,0, by=.5))#cover the mean line with a grey linelines(x=object[,1]-max(object[,1]),y=object[,3], type="l", lwd=2, col="grey60")#plot empirical LTTobj2<-ltt(tree, plot=FALSE)#flip x axis values to negativeobj2$times<-obj2$times-max(obj2$times)#plot on top of existing plotplot(obj2, add=TRUE, log.lineages=FALSE, col="black", lwd=4, lty=1)```## Check whether early burst is found in each replicate: Rep1(r3)```{r}#read in the consensus treetree <- ape::read.nexus(file ="~/Desktop/nmel.ceyx.rad/snapp/rep1.r3/rep1.con.tree")#remove root so that only the diversification rate within the northern Melanesian clade is assessedtree<-drop.tip(tree, "margarethae")#set h = the height of the empirical treeh<-max(nodeHeights(tree))h#set x as 1000 evenly spaced numbers between zero and the height of the input treex<-seq(0,h,by=h/1000)#set birth rate equal to a log-linear rate, following a pure birth predictionb<-(log(Ntip(tree))-log(2))/h#simulate 1000 trees with the same height and tip number as our empirical tree, plus a log-linear diversification ratetrees<-pbtree(b=b, n=Ntip(tree), t=h,nsim=1000 ,method="direct", quiet=TRUE)#calculate 95% CI for the 1000 trees simmed above under a pure birth/death modelobject<-ltt95(trees,log=TRUE,plot=FALSE, mode ="mean")#set up the plotpar(bty="l")#plot background 95% CI based on simulated dataplot(object, labels=FALSE, shaded=TRUE, xaxis ="negative", bg="grey90")axis(side =1, at=seq(-4,0, by=.5))#cover the mean line with a grey linelines(x=object[,1]-max(object[,1]),y=object[,3], type="l", lwd=2, col="grey60")#plot empirical LTTobj2<-ltt(tree, plot=FALSE)#flip x axis values to negativeobj2$times<-obj2$times-max(obj2$times)#plot on top of existing plotplot(obj2, add=TRUE, log.lineages=FALSE, col="black", lwd=4, lty=1)```## Check whether early burst is found in each replicate: Rep2(r1)```{r}#read in the consensus treetree <- ape::read.nexus(file ="~/Desktop/nmel.ceyx.rad/snapp/rep2.con.tree")#remove root so that only the diversification rate within the northern Melanesian clade is assessedtree<-drop.tip(tree, "margarethae")#set h = the height of the empirical treeh<-max(nodeHeights(tree))h#set x as 1000 evenly spaced numbers between zero and the height of the input treex<-seq(0,h,by=h/1000)#set birth rate equal to a log-linear rate, following a pure birth predictionb<-(log(Ntip(tree))-log(2))/h#simulate 1000 trees with the same height and tip number as our empirical tree, plus a log-linear diversification ratetrees<-pbtree(b=b, n=Ntip(tree), t=h,nsim=1000 ,method="direct", quiet=TRUE)#calculate 95% CI for the 1000 trees simmed above under a pure birth/death modelobject<-ltt95(trees,log=TRUE,plot=FALSE, mode ="mean")#set up the plotpar(bty="l")#plot background 95% CI based on simulated dataplot(object, labels=FALSE, shaded=TRUE, xaxis ="negative", bg="grey90")axis(side =1, at=seq(-4,0, by=.5))#cover the mean line with a grey linelines(x=object[,1]-max(object[,1]),y=object[,3], type="l", lwd=2, col="grey60")#plot empirical LTTobj2<-ltt(tree, plot=FALSE)#flip x axis values to negativeobj2$times<-obj2$times-max(obj2$times)#plot on top of existing plotplot(obj2, add=TRUE, log.lineages=FALSE, col="black", lwd=4, lty=1)```## Check whether early burst is found in each replicate: Rep2(r2)```{r}#read in the consensus treetree <- ape::read.nexus(file ="~/Desktop/nmel.ceyx.rad/snapp/rep2.r2/rep2.con.tree")#remove root so that only the diversification rate within the northern Melanesian clade is assessedtree<-drop.tip(tree, "margarethae")#set h = the height of the empirical treeh<-max(nodeHeights(tree))h#set x as 1000 evenly spaced numbers between zero and the height of the input treex<-seq(0,h,by=h/1000)#set birth rate equal to a log-linear rate, following a pure birth predictionb<-(log(Ntip(tree))-log(2))/h#simulate 1000 trees with the same height and tip number as our empirical tree, plus a log-linear diversification ratetrees<-pbtree(b=b, n=Ntip(tree), t=h,nsim=1000 ,method="direct", quiet=TRUE)#calculate 95% CI for the 1000 trees simmed above under a pure birth/death modelobject<-ltt95(trees,log=TRUE,plot=FALSE, mode ="mean")#set up the plotpar(bty="l")#plot background 95% CI based on simulated dataplot(object, labels=FALSE, shaded=TRUE, xaxis ="negative", bg="grey90")axis(side =1, at=seq(-4,0, by=.5))#cover the mean line with a grey linelines(x=object[,1]-max(object[,1]),y=object[,3], type="l", lwd=2, col="grey60")#plot empirical LTTobj2<-ltt(tree, plot=FALSE)#flip x axis values to negativeobj2$times<-obj2$times-max(obj2$times)#plot on top of existing plotplot(obj2, add=TRUE, log.lineages=FALSE, col="black", lwd=4, lty=1)```## Check whether early burst is found in each replicate: Rep2(r3)```{r}#read in the consensus treetree <- ape::read.nexus(file ="~/Desktop/nmel.ceyx.rad/snapp/rep2.r3/rep2.con.tree")#remove root so that only the diversification rate within the northern Melanesian clade is assessedtree<-drop.tip(tree, "margarethae")#set h = the height of the empirical treeh<-max(nodeHeights(tree))h#set x as 1000 evenly spaced numbers between zero and the height of the input treex<-seq(0,h,by=h/1000)#set birth rate equal to a log-linear rate, following a pure birth predictionb<-(log(Ntip(tree))-log(2))/h#simulate 1000 trees with the same height and tip number as our empirical tree, plus a log-linear diversification ratetrees<-pbtree(b=b, n=Ntip(tree), t=h,nsim=1000 ,method="direct", quiet=TRUE)#calculate 95% CI for the 1000 trees simmed above under a pure birth/death modelobject<-ltt95(trees,log=TRUE,plot=FALSE, mode ="mean")#set up the plotpar(bty="l")#plot background 95% CI based on simulated dataplot(object, labels=FALSE, shaded=TRUE, xaxis ="negative", bg="grey90")axis(side =1, at=seq(-4,0, by=.5))#cover the mean line with a grey linelines(x=object[,1]-max(object[,1]),y=object[,3], type="l", lwd=2, col="grey60")#plot empirical LTTobj2<-ltt(tree, plot=FALSE)#flip x axis values to negativeobj2$times<-obj2$times-max(obj2$times)#plot on top of existing plotplot(obj2, add=TRUE, log.lineages=FALSE, col="black", lwd=4, lty=1)```## Check whether early burst is found in each replicate: Rep3(r1)```{r}#read in the consensus treetree <- ape::read.nexus(file ="~/Desktop/nmel.ceyx.rad/snapp/rep3.consensus.tree")#remove root so that only the diversification rate within the northern Melanesian clade is assessedtree<-drop.tip(tree, "margarethae")#set h = the height of the empirical treeh<-max(nodeHeights(tree))h#set x as 1000 evenly spaced numbers between zero and the height of the input treex<-seq(0,h,by=h/1000)#set birth rate equal to a log-linear rate, following a pure birth predictionb<-(log(Ntip(tree))-log(2))/h#simulate 1000 trees with the same height and tip number as our empirical tree, plus a log-linear diversification ratetrees<-pbtree(b=b, n=Ntip(tree), t=h,nsim=1000 ,method="direct", quiet=TRUE)#calculate 95% CI for the 1000 trees simmed above under a pure birth/death modelobject<-ltt95(trees,log=TRUE,plot=FALSE, mode ="mean")#set up the plotpar(bty="l")#plot background 95% CI based on simulated dataplot(object, labels=FALSE, shaded=TRUE, xaxis ="negative", bg="grey90")axis(side =1, at=seq(-4,0, by=.5))#cover the mean line with a grey linelines(x=object[,1]-max(object[,1]),y=object[,3], type="l", lwd=2, col="grey60")#plot empirical LTTobj2<-ltt(tree, plot=FALSE)#flip x axis values to negativeobj2$times<-obj2$times-max(obj2$times)#plot on top of existing plotplot(obj2, add=TRUE, log.lineages=FALSE, col="black", lwd=4, lty=1)```## Check whether early burst is found in each replicate: Rep3(r2)```{r}#read in the consensus treetree <- ape::read.nexus(file ="~/Desktop/nmel.ceyx.rad/snapp/rep3.r2/rep3.con.tree")#remove root so that only the diversification rate within the northern Melanesian clade is assessedtree<-drop.tip(tree, "margarethae")#set h = the height of the empirical treeh<-max(nodeHeights(tree))h#set x as 1000 evenly spaced numbers between zero and the height of the input treex<-seq(0,h,by=h/1000)#set birth rate equal to a log-linear rate, following a pure birth predictionb<-(log(Ntip(tree))-log(2))/h#simulate 1000 trees with the same height and tip number as our empirical tree, plus a log-linear diversification ratetrees<-pbtree(b=b, n=Ntip(tree), t=h,nsim=1000 ,method="direct", quiet=TRUE)#calculate 95% CI for the 1000 trees simmed above under a pure birth/death modelobject<-ltt95(trees,log=TRUE,plot=FALSE, mode ="mean")#set up the plotpar(bty="l")#plot background 95% CI based on simulated dataplot(object, labels=FALSE, shaded=TRUE, xaxis ="negative", bg="grey90")axis(side =1, at=seq(-4,0, by=.5))#cover the mean line with a grey linelines(x=object[,1]-max(object[,1]),y=object[,3], type="l", lwd=2, col="grey60")#plot empirical LTTobj2<-ltt(tree, plot=FALSE)#flip x axis values to negativeobj2$times<-obj2$times-max(obj2$times)#plot on top of existing plotplot(obj2, add=TRUE, log.lineages=FALSE, col="black", lwd=4, lty=1)```## Check whether early burst is found in each replicate: Rep3(r3)```{r}#read in the consensus treetree <- ape::read.nexus(file ="~/Desktop/nmel.ceyx.rad/snapp/rep3.r3/rep3.con.tree")#remove root so that only the diversification rate within the northern Melanesian clade is assessedtree<-drop.tip(tree, "margarethae")#set h = the height of the empirical treeh<-max(nodeHeights(tree))h#set x as 1000 evenly spaced numbers between zero and the height of the input treex<-seq(0,h,by=h/1000)#set birth rate equal to a log-linear rate, following a pure birth predictionb<-(log(Ntip(tree))-log(2))/h#simulate 1000 trees with the same height and tip number as our empirical tree, plus a log-linear diversification ratetrees<-pbtree(b=b, n=Ntip(tree), t=h,nsim=1000 ,method="direct", quiet=TRUE)#calculate 95% CI for the 1000 trees simmed above under a pure birth/death modelobject<-ltt95(trees,log=TRUE,plot=FALSE, mode ="mean")#set up the plotpar(bty="l")#plot background 95% CI based on simulated dataplot(object, labels=FALSE, shaded=TRUE, xaxis ="negative", bg="grey90")axis(side =1, at=seq(-4,0, by=.5))#cover the mean line with a grey linelines(x=object[,1]-max(object[,1]),y=object[,3], type="l", lwd=2, col="grey60")#plot empirical LTTobj2<-ltt(tree, plot=FALSE)#flip x axis values to negativeobj2$times<-obj2$times-max(obj2$times)#plot on top of existing plotplot(obj2, add=TRUE, log.lineages=FALSE, col="black", lwd=4, lty=1)```